# -*- coding: utf-8 -*-
##DESCRIPTION: Analysis of landslide runout 
##
##References:
##
##[1] Hungr, O., Corominas, J. & Eberhardt, E. (2005). Estimating landslide motion mechanism, travel distance and velocity.
##In Landslide Risk Management. Edited by O. Hungr, R. Fell, R. Couture & E. Eberhardt, A.A. Balkema, Leiden, 99-128.

##[2] Corominas, J. 1996. The angle of reach as a mobility index for small and large landslides.
##Canadian Geotechnical Journal33: 260-271.

##[3] Hunter, G. & Fell, R. 2003. Travel distance angle for rapid-landslides in constructed and natural soil slopes.
##CanadianGeotechnical Journal 40(6): 1123-1141.

##[4] Hurlimann, M., Medina, V., Bateman, A., Copons, R., and Altimir, J. 2007. Comparison of different techniques to analyse
##the mobility of debris flows during hazard assessment—Case study in La Comella catchment, Andorra.
##In Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, Chen & Major, eds. Millpress, Netherlands.



#modules and library
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
from scipy.interpolate import UnivariateSpline


#functions
def compute_HoL(V,A,B):
    #log (H/L) = A +  B* log(V)    [1] eq. 2 and 3 
    #V in m3
    HoL=10**(A+B*np.log10(V))
    return -HoL

def beta_slope_profile(numprofiles, max_a, min_a, max_b, min_b, elev,base):
    #generating random parameters for beta distribution
    a_list=np.random.rand(numprofiles)*(max_a-min_a)+min_a
    b_list=np.random.rand(numprofiles)*(max_b-min_b)+min_b
    b_list=a_list/b_list
    x=np.linspace(0,base,numpts)
    y_list=[]
    for j in range(len(a_list)):
        #3. Generating topo profiles
        y=beta.cdf(x,a_list[j],b_list[j],loc=0,scale=base)
        y_list.append(y[::-1])
    return x,y_list

#inputs
#Table 4 in [1] from [2]
LStype=['Rockfalls', 'Translational slides', 'Debris flows', 'Earth flows']
Paths=['All', 'Obstructed', 'Unobstructed']
A=[0.21,0.231,0.167,
   -0.159,-0.133,-0.143,
   -0.012,-0.031,-0.049,
   -0.214,np.nan,-0.220]
B=[-0.109,-0.091,-0.119,
   -0.068,-0.057,-0.080,
   -0.105,-0.102,-0.108,
   -0.070,np.nan,-0.138]
R2=[0.76,0.83,0.92,
    0.67,0.76,0.80,
    0.76,0.87,0.85,
    0.65,np.nan,0.91]
indx=pd.MultiIndex.from_product((LStype,Paths), names=['Landslide type','Path'])
tab4=pd.DataFrame(data={'A':A,'B':B,'R2':R2}, index=indx)
tab4=tab4.reset_index()
print tab4



#1. plotting V vs H/L
plt.figure(1)
V=np.logspace(2,6, 5*(6-2)+1)  #m3
LStype_list=[2,5,8,11]

for i in LStype_list:
    LStype=tab4['Landslide type'][i]+","+tab4['Path'][i]
    A=tab4['A'][i]
    B=tab4['B'][i]
    HoL=compute_HoL(V,A,B)
    plt.plot(V,-HoL,'.-', label=LStype)
plt.legend()
plt.semilogx()
plt.xlabel('Volume, m^3')
plt.ylabel('H/L')
#plt.show()

#2. generating topographic profiles
#2.1 define number of profiles
numprofiles=2000
#2.2 defining beta distribution parameters for topo profile
max_a=7
min_a=2
max_b=5  #increase max_b,min_b to produce steeper slopes at the source, as desired
min_b=2
#2.3 defining basic topo profile
elev=435   #height of source above base
base=1360  #length of base
numpts=100 #number of points along profile to generate
#2.4 generating profiles
x,y_list=beta_slope_profile(numprofiles, max_a, min_a, max_b, min_b, elev,base)


#3.2 defining values of volume to analyze and plot
V=np.asarray([5,10,50,100])*10**3 #m^3

#3.3 defining landslide type to analyze
for h in LStype_list:
    #3.4 selecting parameters for current landslide type
    LStype=tab4['Landslide type'][h]
    A=tab4['A'][h]
    B=tab4['B'][h]
    print tab4[(tab4.index==h)]

    #3.5 creating plot and subplots for current landslide type and different volumes
    fig,ax=plt.subplots(nrows=len(V), ncols=1, sharex=True, sharey=True)
    fig.suptitle(LStype)
    for i in range(len(V)):
        plt.sca(ax[i])

        #3.6 plot basic planar profile
        plt.plot([0,base],[elev,0],'k--',lw=2)

        #3.7 plotting random topographic profiles
        v=V[i]
        HoL=compute_HoL(v,A,B)
        for j in range(numprofiles):
            y=y_list[j]
            s=UnivariateSpline(x,y*elev,s=0)
            #3.8 plotting admissible random profiles
            if s(0,1)<HoL:
                plt.plot(x,s(x),'k-',lw=0.2,color='0.2')

        #3.9 plotting H/L profile
        H=HoL*x+elev
        plt.plot([0,(-elev/HoL)],[elev,0], 'r-', lw=2,label="H/L="+str(round(-HoL,2))+"; V="+str(v)+" m^3")

        #3.10 formatting axes
        plt.axis('equal')
        plt.legend(fontsize='x-small')
        if v==max(V):
            plt.xlabel("horizontal distance, m")
            
plt.show()
        
    
    
               




###Table 5 in [1] from [3]
##Paths=['Unconfined', 'Partly confined', 'Confined']
##A=[0.77,  0.69, 0.54]
##B=[0.087, 0.110,0.27]
##R2=[0.71, 0.52,0.85]
##SD=[0.095,0.11,0.027]
##tab5=pd.DataFrame(data={'A':A, 'B':B, 'R2':R2, 'SD':SD}, index=[Paths])
##
##tab4['V,10^3 m^3']=V*np.ones(len(tab4))
##tab4['H, m']=H*np.ones(len(tab4))
##tab4['L, m']=compute_L(tab4['H, m'],tab4['V,10^3 m^3'],tab4['A'],tab4['B'])
###tab4['H/tan32']=tab4['H']/np.tan(np.radians(32))
###tab4['LE']=tab4['L']-tab4['H/tan32']
##print tab4
##
##tab5['V']=V*np.ones(len(tab5))
##tab5['H']=H*np.ones(len(tab5))
##tab5['L']=compute_L(tab5['H'],tab5['V'],tab5['A'],tab5['B'])
##
##
